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INTRODUCTION  AND  SUMMARY 


Multitarget;  Tracking  Studies  (MTS)  is  a  research  effort  which  has 
the  objectives  of  developing  and  evaluating  a  new  concept  for  tracking 
multiple  targets.  The  algorithms  developed  in  this  program  (which  will 
be  referred  to  as  the  MTS  algorithm  or  processor)  will  complement  and 
enhance  currently  used  tracking  techniques.  While  the  main  goal  of  the 
program  is  directed  towards  multiple  targets,  MTS  is  expected  to  have  a 
significant  impact  on  the  single  target  case  as  well. 

In  this  report,  we  summarize  the  main  theoretical  developments  and 
some  preliminary  performance  evaluation  results.  This  work  is  part  of 
the  first  phase  of  the  MTS  research  program  in  which  the  single  target 
case  was  studied.  Results  so  far  have  been  very  promising  and  we  expect 
to  adapt  the  techniques  developed  in  this  initial  phase  to  handling 
multiple  targets  during  the  second  phase  of  the  program. 

An  important  point  here  is  the  following.  It  has  not  been  the 
objective  of  this  study  to  surpass  conventional  estimation  performance 
of  spectrum  and  TDOA  analyzers.  The  original  emphasis  was  upon  demon¬ 
stration  that  a  framework  in  which  these  functions  can  be  carried  out 
naturally  for  multiple  targets  is  one  in  which  performance  on  individual 
targets  can  be  maintained.  It  was  sufficient,  therefore,  to  demonstrate 
that  even  for  low  SNR  (-10dB-0dB)  the  MTS  performance  compares  with 
conventional  approaches.  These  approaches  make  their  processing  gains 
by  integration,  a  procedure  also  available  to  MTS.  Because  pre- integration 
performance  of  MTS  was  so  encouraging,  no  further  pursuit  toward  com¬ 
parison  was  undertaken.  However,  it  turned  out  that  the  new  approach 
to  adaptive  signal  processing  implicit  in  the  MTS  processing  could  be 
developed  to  one  offering  substantial  improvement  over  current  techniques. 
While  our  research  effort  is  directed  towards  the  development  of 
tracking  algorithms,  the  signal  modeling  approach  has  a  much  wider 
applicability  to  the  Navy's  signal  processing  problems.  In  particular, 
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the  MTS  algorithm  is  directly  applicable  to  a  number  of  adaptive 
signal  processing  problems,  including:  line  enhancement,  high  resolution 
spectral  estimation,  noise  cancelling  and  channel  equalization.  Applying 
the  proposed  ARMA  modeling  techniques  to  some  of  these  problems  has 
already  resulted  in  substantial  performance  improvements.  The  analysis 
of  the  MTS  algorithm  as  it  is  used  for  adaptive  signal  processing  is 
summarized  in  [10]. 

The  MTS  concept  is  based  on  modeling  the  observed  data  as  an  ARMA 
process.  The  parameters  of  the  model  provide  a  compact  representation 
of  target  parameters  such  as  spectrum  and  TDOA/bearing.  These  parameters 
can,  therefore,  be  used  as  inputs  to  a  tracking  algorithm,  a  target 
classification  program,  etc.  Section  2  of  the  report  describes  the  MTS 
concept  and  how  it  fits  into  an  overall  system. 

A  major  part  of  our  effort  has  been  directed  towards  developing  and 
coding  the  basic  MTS  algorithm.  This  algorithm  is  a  parameter  estimation 
technique  which  recursively  computes  a  set  of  ARMA  parameters  from  the 
observed  data  sequence.  This  algorithm  is  now  implemented  as  an  inter¬ 
active  computer  program  and  it  provides  a  very  powerful  and  flexible  signal 
processing  tool.  This  program  will  be  the  core  of  our  future  MTS  work. 
Section  3  of  the  report  describes  the  algorithm  and  its  main  features. 

The  main  issues  addressed  so  far  are  TDOA  estimation  and  estimation  of 
the  spectral  parameters  of  the  target  under  different  signal- to- noise  ratio 
conditions.  Several  synthetic  test  cases,  both  narrowband  and  broadband, 
were  used  to  evaluate  the  performance  of  the  MTS  algorithm.  Results  were 
very  encouraging. 

For  high  SNR  (20dB  and  above)  the  algorithm  provided  excellent  results, 
and  had  no  problems  in  converging  to  the  right  spectral/TDOA  parameters. 

In  moderate  SNR  (0-20dB) ,  serious  convergence  problems  were  initially 
experienced. 
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A  significant  amount  of  effort  was  devoted  to  studying  and  solving 
these  problems.  Our  solution  provides  an  important  step  in  extending  the 
range  of  applicability  of  recursive  parameter  estimation  algorithms  to  low 
SNR  situations.  A  number  of  publications  on  this  topic  are  in  preparation. 
Currently,  we  are  able  to  get  good  spectral  and  TDOA  estimates  for  SNR's  in 
the  0-20dB  range.  Our  experience  with  low  SNR  (-10dB-0dB)  has  been  that 
performance  matching  or  exceeding  conventional  approaches  can  be  achieved. 
No  special  difficulties  were  observed  at  low  SNR,  but  we  feel  that  some 
refinements  of  the  algorithm  may  improve  performance  even  further. 

The  positive  results  obtained  so  far  will  provide  the  basis  for  our 
next  phase  of  research,  in  which  the  multitarget  tracking  algorithm  will 
be  developed.  Our  research  will  be  performed  in  two  steps: 

(i)  Complete  the  single  target  tracking  algorithm. 

Here,  we  will  concentrate  on  the  issues  associated  with  the  operation 
of  the  MTS  algorithm  at  low  SNR.  We  will  make  the  algorithm  more  robust 
by  pre-filtering  and  other  methods,  test  its  tracking  capability  on 
synthetic  data  with  time  varying  parameters  and  develop  performance 
bounds  to  evaluate  its  performance  against  suitable  standards. 

(ii)  Development  and  testing  of  multitarget  algorithms  for  the 
high  SNR  case. 

Here  we  will  develop  and  evaluate  a  candidate  algorithm  for  tracking 
several  targets.  The  objective  will  be  to  demonstrate  the  capability 
of  an  MTS  algorithm  to  provide  consistent  tracks  for  several  targets. 

In  particular,  we  will  investigate  the  special  structural  properties 
of  multi-input  multi-output  (MIMO)  systems  of  the  type  used  to  model  the 
multitarget  tracking  problem,  and  study  questions  of  identif iability 
and  uniqueness.  The  extension  of  the  MTS  approach  to  tracking  multiple 
targets  at  low  SNR  will  be  deferred  to  a  third  phase  of  the  project. 
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Several  issues  need  to  be  investigated  in  order  to  achieve  such  an 
extension,  including:  the  convergence  of  the  MIMO  RML  algorithm, 
development  of  pre-filtering  and  other  mechanisms  for  improved  con¬ 
vergence  and  analysis  of  the  uniqueness  and  identif iability  issues 
under  low  SNR  conditions. 

This  plan  of  work  is  summarized  in  the  following  schematic: 


SINGLE  MULTIPLE 

TARGET  TARGETS 


PHASE  I 

PHASE  II 

High  SNR 

— 

Low  SNR 

(i) 

f 

_ _ 

f 

PHASE  II 

_ EHASE-LLI 

We  believe  that  the  next  phase  of  the  MTS  project  will  result  in 
significant  contributions  to  the  areas  of  multitarget  tracking,  adaptive 
signal  processing,  multichannel  parameter  estimation  and  modeling  of 
vector  time-series. 
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2.  SYSTEM  DESCRIPTION 

The  MTS  algorithm  is  a  coherent,  time-domain  signal  processing  tech¬ 
nique  for  extracting  target  parameters  (spectrum  and  TDCA/bearing)  from 
multisensor  data.  The  sensors  may  be  the  elements  of  one  or  several  arrays. 
The  MTS  algorithm  may  operate  directly  on  wideband  sensor  data,  as  depicted 
in  Figure  1. 


TARGET 
PARADE i t ao 


However,  since  the  computational  requirements  of  the  MTS  algorithm 
increase  in  proportion  to  the  bandwidth  of  the  input  signal  (as  will  be  show 
later),  it  is  desirable  to  reduce  the  bandwidth  of  the  sensor  signals  before 
handing  them  to  the  MTS  algorithm.  This  can  be  done  by  a  preprocessing  step, 
in  which  one  or  more  spectral  bands  of  interest  are  shifted  in  frequency  to 
provide  a  combined,  relatively  narrowband  signal,  as  depicted  in  Figure  2. 


A 


0  B 


Figure  2:  Bandwidth  Selection 
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This  bandwidth  reduction  can  be  implemented  in  many  different  ways. 
A  block  diagram  of  one  possible  implementation  is  depicted  in  Figure  3. 
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a)  A  Single  Spectral  Band  (Centered  at  f  ) 
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t 
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FILTER  k 
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TO 

MTS 

PROCESSOR 


b)  Multiple  Spectral  Bands 


Figure  3.  Preprocessing  for  Bandwidth  Reduction 


In  the  rest  of  this  report,  we  will  always  assume  that  this  prepro¬ 
cessing  step  has  been  performed,  and  that  the  MTS  processor  is  handed  data 
with  a  bandwidth  of  B  Hz.  The  data  is  sampled  at  the  Nyquist  rate,  thus 

AT  *  sampling  interval  =  1/2B  (1) 


A  typical  spectrum  of  the  signal  at  the  input  of  the  MIS  algorithm 
will  contain  several  spectral  lines  in  a  noise  background,  as  depicted  in 
Figure  4.  This  spectrum  was  obtained  by  performing  an  N  point  FFT  of 


the  data  where 


2  SINE  WAVES  IN  NOISE  —  SNR=0  dB 


N  =  512 

Figure  4:  Typical  Spectrum  of  Synthetic  Data 


N  =  number  of  data  points, 
which  corresponds  to  an  integration  time  of 
T  =  NAT  =  N/2B. 


The  frequency  resolution  of  this  spectral  plot  is  Af  Hz  per  point 


Thus,  the  i-th  point  on  the  plot  represents  a  frequency  f^,  where 


f 


i 


i  _  21B 
NAT  "  N  ' 


(4) 


In  this  report,  we  will  define  the  signal- to-noise  ratio  (SNR)  as 
the  ratio  of  the  total  signal  energy  to  the  total  noise  energy  in  the  band¬ 
width  B  which  is  provided  to  the  MTS  processor.  The  signal  and  noise 
processes  are  generated  by  a  synthetic  data  generator  which  produces  for 
each  sensor  a  data  sequence  y^t). 

yf(t)  =  si(t)  +  nt(t)  (5) 

si(t)  ■  signal  arriving  at  sensor  i 

n^(t)  =  measurement  noise  at  sensor  i  (white  Gaussian 
noise,  independent  from  sensor  to  sensor) . 


The  total  signal  and  noise  energies 


(si*  Ni^  are  comPuted 

(6a) 

(6b) 


and  the  corresponding  signal-to-noise  ratio  is  given  by 

SNRi  =  Si/Ni  (7) 

The  noise  energy  is  related  to  its  spectral  power  density  by 

N  =■  N  B  (8) 

o 

where  N  is  noise  energy  per  Hz. 
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The  MTS  algorithm  is  based  on  the  idea  of  fitting  an  autoregres¬ 
sive  moving-average  (ARMA)  model  to  the  observed  time  series  (see 
Appendix  A  for  a  more  detailed  explanation).  The  basic  model  is  depicted 
in  Figure  5.  The  autoregressive  (AR)  part  of  the  model  provides  informa¬ 
tion  about  the  spectrum  of  the  target,  while  the  moving-average  (MA)  part 
gives  the  TDOA  information.  Thus,  once  an  ARMA  model  has  been  fit  to  the 
observed  data,  all  the  target  parameters  can  be  obtained  from  the  ARMA 
coefficients  (a^.B^}.  The  spectral  estimate  of  the  target  can  be 
obtained  by  an  FFT  of  the  impulse  response  of  the  AR  model  portion*. 

More  precisely,  we  can  FFT  the  time  series  x(t) 
n 

a 

x(t)  =  -  £  +  u(t)  (9) 

i=l 


u(t) 


t=0 

t#0 


The  normalized  estimate  of  the  target  power  spectrum  Sx(i)  is 
given  by 


S  (i)  = 
x 


(10) 


where  {x^}  is  the  FFT  of  {x(t)}.  The  interpretation  of  the  frequency 
corresponding  to  the  i-th  spectral  estimate  (S^ (i) )  is  given  by  (4). 

The  TDOA  estimates  can  be  obtained  by  looking  at  the  b^  coeffi¬ 
cients,  as  depicted  in  Figure  6.  The  TDOA  corresponding  to  a  difference 
of  one  (in  order)  is  AT,  where  AT  is  the  sampling  rate  of  the  data  at 
the  input  to  the  MTS  algorithm.  This  is  not  necessarily  the  ultimate 
resolution  of  our  TDOA  estimation,  since  finer  resolution  can  be  achieved 
by  interpolation.  A  more  detailed  discussion  of  this  point  can  be  found 
in  Appendix  B. 


*Note  that  we  could  also  evaluate  a  z-transform. 


c.j  =  Attenuation, 


D-  =  Delay 


Target: 

x(t)  =  -  5^  ai  x(t-i)  +  u(t) 


X  ( z )  = 


U(z) 


Receiver: 


y(t)  = 


c-|X(t-0-| ) 
c2x(t-D2) 

Y(z)  = 

•°1 

c,z 

-D 

u2 

C2Z 

f3x(t“D3)_ 

.c32"°i 

B(z) 


The  Overall  Model 
n 


m 


y(t)  =  -  22  a.y(t-i)  +  22  B.u(t-i) 

i=l  1  i=l  1 


V(z)  =  III}  U(z) 


Figure  5.  An  ARMA  Model  for  the  Single  Target  Case 


X(z) 
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Figure  6.  Estimating  TDOA-s  From  the  tK  Coefficients 

A  change  of  AT  in  the  TDOA  can  be  translated  into  a  change  of 
bearing  A0  by  (see  Figure  7) , 

AS  =  vAT/Lcos0  *  v/2BLcos9 

where 

L  =  distance  between  the  two  sensors 
9  *  bearing 
v  »  sound  velocity. 


SENSOR  1 


TDOA  =  L  sine/v 
AT  =  L  COS0  A0/v 

Figure  7:  The  Geometry  for  TDOA  Computation 


The  total  angular  extent  over  which  the  MTS  processor  can  be  "steered" 
is  given  by  where  n^  =  the  number  of  b  coefficients,  i.e., 

the  order  of  the  MA  model.  This  angular  extent  can,  of  course,  be  increased 
by  removing  bulk  delays  prior  to  the  MTS  processor.  If  necessary,  several 
MTS  processors  can  be  run  in  parallel,  each  covering  a  section  of  Ac-n^ 
degrees. 

Consider,  for  example,  the  following  representative  case: 

L  =  1200  meters 
v  =  1490  m/sec 

%  “  20 
B  -  10  Hz 

then:  A0  =  3.55° 

A0*n^  *  71°. 

The  presence  of  doppler  shifts  in  the  received  signals  will  be  handled 
in  the  MTS  algorithm  by  computing  a  different  set  of  { a^}  for  each  sensor. 
In  other  words,  different  sensors  will  observe  different  (shifted)  spectral 
lines.  This  feature  of  the  algorithm  has  not  been  tested  yet,  but  more 
details  can  be  found  in  Section  3. 


p 


The  estimated  target  parameters  computed  by  the  algorithm  will  be 
used  as  an  input  to  various  post  processors  which  extract  operational 
parameters  such  as  target  location  (coordinates),  target  signature,  target 
type,  etc.  An  overall  block  diagram  of  the  processing  in  an  MTS  system  is 
depicted  in  Figure  8. 


Post  Process1n9 


Figure  8.  Block  Diagram  of  a  Basic  MTS  System 
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3. 


THE  MTS  ALGORITHM 


The  core  of  Che  MTS  processor  is  a  recursive  parameter  estimation 
algorithm  which  estimates  ARMA  coefficients  from  an  observed  data  sequence. 
Algorithms  of  this  type  have  been  developed  in  the  context  of  adaptive  con¬ 
trol  [1].  Our  application  of  this  class  of  algorithms  to  acoustic  signal 
processing  seems  to  be  a  pioneering  effort  which  promises  to  lead  to  a 
whole  new  class  of  adaptive  signal  processing  techniques.  Some  important 
modifications  are  required  in  transforming  this  type  of  algorithm  from  the 
control  context  to  the  signal  processing  context.  A  significant  part  of 
our  research  effort  was  directed  to  investigation  and  development  of  these 
modifications.  A  key  development,  which  is  described  later  in  this  section, 
was  the  Improvement  of  the  convergence  properties  of  the  algorithm. 

The  Basic  Algorithm 

Several  versions  of  recursive  parameter  estimation  algorithms  have 
been  coded  and  tested: 

(1)  Recursive  Least  Squares  (RLS) 

(2)  Recursive  Maximum  Likelihood  (RML1) 

(3)  Modified  Recursive  Maximum  Likelihood  (RMLP) 

(4)  Recursive  Maximum  Likelihood  with  Prefiltering  (RML2) 

Initial  experiments  indicated  that  the  RML2  algorithm  is  most  suitable 
for  our  application.  We  will  therefore  describe  here  only  the  RML2  algorithm. 
For  a  more  detailed  description  of  all  of  these  algorithms  see  [2],  [3]. 

The  RML2  algorithm  estimates  the  parameters  of  an  ARMA  model  of  the 
following  type: 

na  nb  *c 

y(t)  *  —  y*  a  y (t)  +  b  u(t-i)  +  V  c.  e(t-i)  (12) 

i=l  1  i-1  1  i-0 

where  e(t)  is  an  (unobservable)  white  noise  process.  The  presence  of  the 
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"*“»  buuuhot  hujj> 


coefficients  enables  us  to  handle  correlated  measurement  noise  and 
the  case  of  unknown  inputs.  It  is  assumed  that  cq=1.  Equation  (12) 
can  be  written  more  compactly  as 

y (t)  *  0T(t)0  +  e(t)  (13) 

where 

4>T(t)  *  [-y(t-l) . -y(t-na);  u(t-l) , . . .  ^(t-i^) ;  e(t-l) , . . . ,e(t-nc) ] 

T 

9  *  [a- ,  *  *  - »  a  ;  b. ,  b  ;  c1 ,  • • • >  c  ] 

I  n  1  n  i  n 

a  a  c 

The  dimension  of  0  and  <J>  is 

n  *  n  +  a.  +  n  . 
a  d  c 

Since  Eq.  (13)  is  linear  in  the  unknown  parameters  (the  components 
of  6),  a  recursive  estimation  algorithm  is  obtained  by  the  following 
set  of  Kalman  filter  equations: 

0(t+l)  =  0(t)  +  K (t+1)  e(t+l)  '(14a; 

K(t+1)  -  P(t)  <Kt+l)/(A  +  4>T(t+l)P(t)4>(t+l))  *  P(t+1)  <Mt+l)  (14b) 

P(t+1)  *  [P(t)-P(t)0(t+l)<f>^(t+l)P(t)/  (X-HJ>T(t+l)P(t)$(t+l))  ]/A  =  (14c) 

■  error  covariance  of  the  parameters. 
e(t+l)  ■  y(t+l)  -  4> ( t+1) T  9(t)  *  prediction  error,  (14d) 

with  initial  conditions, 

P(o)  «  ol,  a  ■  a  scalar  parameter 

A. 

0(o)  *  0  or  8  ,  a  prior  estimate. 
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The  parameter  X  represents  data  windowing,  i.e.,  it  is  the  "forgetting 
factor"  of  the  algorithm.  Various  (time-varying  as  well  as  fixed)  values 
of  this  parameter  have  been  tried  out.  To  facilitate  the  convergence  of 
the  algorithm  on  short  data  sequences  the  following  was  found  to  give  the 
best  results: 

X (t+1)  =  X  X(t)  +  (1-X  )  (15) 

o  o 

where  X(o),  Xq  are  specified  parameters. 

Other  choices  for  X  are  described  in  [2].  Two  additional  quantities 
which  are  useful  to  keep  track  of  the  numerical  behavior  of  the  algorithm 
are: 

n 

trace  (P (t)  }  =  £  P,.  (t)  (16) 

i=l  1 

and 

r|(t)  =  — y-  trace  (P ( t ) >  trace  (P(t)-1}  =  (17) 

n 

*  a  measure  of  how  close  to  singular  is  P(t). 

The  only  difficulty  with  the  algorithm  described  above  is  that  it 
requires  knowledge  of  the  unobservable  noise  sequence  e(t)  (which  is 
required  for  $(t)).  Since  e(t)  is  unknown,  it  needs  to  be  replaced 
by  some  estimate  of  e(t).  Different  versions  of  the  Recursive  Maximum 
Likehihood  algorithm  are  obtained  by  different  choices  of  the  estimate  of 
e(t).  For  example: 


RML2 


e(t)  -  y(c)  -  MC)T  $(t)  (19) 

In  RML2,  the  unknown  e(t)  is  replaced  by  a  filtered  version  of  the 
prediction  error  e(t).  This  filtering  is  crucial  to  the  proper  conver¬ 
gence  of  the  algorithm  in  MTS  applications. 

The  filtering  is  accomplished  by  replacing  the  He)  vector  which 
is  used  in  Equations  (14c),  (14d)  by  a  version  of  Ht)  filtered  by  1/D(z) 
where 

D(z)  =  l+d,:'^  ...  +  d„  z"nd  (20) 

1  nd 

Sunmary  of  the  Filtering  Equations 

Let  n  =  max  in  ,  n,  ,  n  l 
max  |  a  h  c( 

Define  the  n  x  n  matrix  D 
max  max 


1 


0 


0 

0 


(21) 


Define  the  n  xl  vectors  x.,  x„,  x.  by  the  following  recursions 
max  123 


y(t) 

0 


x^ (t+1)  -  Dv^(t)  + 


0 


x1(o)  -  0 


(22a) 
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x2(t+l)  =  Dx,,(t)  - 


,  x2(o)  -  0 


x3(t+l)  =  Dx3(t)  - 


,  x3(o)  =  0 


Xx(t) 

be 

n  xl 
a 

consisting  of 

the 

first 

n. 

a 

entries 

of 

x^t) 

x2(t) 

be 

V1 

consisting  oi 

the 

first 

% 

entries 

of 

x2(t) 

x3(t) 

be 

n  xl 
c 

consisting  of 

the 

first 

n 

c 

entries 

of 

x3(t) 

*T(t)  «*  [x^(C),  x2(t),  x3(t)] 


The  significance  of  this  filtering  has  to  do  with  the  convergence 
properties  of  recursive  parameter  estimation  algorithms.  Convergence 
analysis  has  shown  ([4]-[6])  that  without  prefiltering,  the  criterion 
for  convergence  is  that 

»<2>  ■  ck  -  i 


be  strictly  positive  real,  i.e., 

Re{H(e-'w)}  >  0  for  all  w. 


Unfortunately,  as  will  be  discussed  later,  this  condition  is  not 
fulfilled  for  general  MTS  signals.  Since  C(z)  is  a  property  of  the 
signal,  and  is  not  under  our  control,  it  is  not  possible  to  guarantee 
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convergence  in  this  case.  Wich  prefiltering,  the  condition  for  conver¬ 
gence  becomes 


H(z) 


D(z) 

C(z) 


strictly  positive  real. 


(25) 


The  choice  of  the  filter  D (z)  is  under  our  control  and,  therefore, 
there  is  hope  of  guaranteeing  convergence.  A  typical  choice  for  D(z) 

[3]  is 


D(z)  =  C(z)  (26) 

The  reasoning  behind  this  choice  is  that  if  C(z)  is  a  good  estimate  of 
C(z),  we  will  get 

“<z)  ‘  -  1 :  1  -  5  =•  °-  <27> 

In  our  preliminary  tests,  we  discovered  that  for  signals  generated 
by  sine  waves  in  white  additive  noise,  this  type  of  filter  was  inadequate 
and  convergence  could  not  be  achieved.  This  problem  was  the  major 
stumbling  block  in  our  initial  research  effort  and  led  to  a  more  careful 
investigation  into  the  convergence  of  RML2  for  MTS  signals.  A  solution 
to  the  problem  has  been  found  and  sucessfully  tested.  The  technique  we 
developed  is  a  significant  contribution  to  the  study  and  application  of 
recursive  parameter  estimation.  The  main  ideas  of  our  technique  are 
described  next. 
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Improved  Pre-Filtering  for  RML2 

To  understand  the  difficulties  inherent  in  the  pre-filtering  problem, 
we  must  first  see  what  the  C(z)  polynomial  means  in  terms  of  the  target 
spectrum  A(z),  the  delay  structure  B(z)  and  the  signal-to-noise  ratio. 

The  observed  signal  y(t)  is  given  by 

y(t)  =  U(t)  +  v(t)  (28) 

signal  measurement 

noise 

i  n 

.  -1  -  a 

1  +  a, z  +...+  a  z 

1  n 

a 

b^  z  +  ...  +  b 

independent  white  noise  processes  with  variance  <j“ 
and  respectively. 

Multiplying  through  by  A(z)  we  get 

A(z)y(t)  =  B(z)u(t)  +  A(z)v(t)  (29) 

Since  neither  u(t)  nor  v(t)  is  directly  measureable,  there  is  no  way 

of  distinguishing  between  them  and  they  can  be  replaced  by  a  white  process 

2 

e(t)  with  variance  such  that 

A(z)y(t)  =  C(z)  e(t) ,  (30) 


where 

A(z)  = 
B(z)  = 
u(t)  ,v(t)  = 


where 

a2c(z)  C(z_1)  =  a2B(z)  B(z-1)  +  a2A(z)  A(z“1).  (31) 

e  u  n 

In  other  words,  C(z)  e(t)  will  have  the  same  spectrum  as  B(z)  u(t)  + 

A(z)  v(t).  To  gain  some  insight  into  what  C(z)  may  look  like, 
we  consider  two  simple  examples.  Both  examples  assume  B(z)  =  b^z  , 
i.e.,  a  pure  delay  propagation  model. 
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(i)  SNR 


a  00 


In  this  case,  a  =0  and  therefore 
n 

2  -12  2 

a  C(z)  c(z  •*•)  =  a 
e  u  D 


or 

C(z)  =  a  b^/a  =  constant  (32) 

u  D  e 

(ii)  SNR  =  0  (or  -oadb) 

In  this  case,  the  second  term  on  the  right-hand  side  of 
dominates,  and  therefore, 

O2  C(z)  C (z"*1)  =  02  A(z)  A(z-1) 
e  n 

or 

C (z)  m  A(z) .  (33) 

In  general,  as  the  SNR  decreases,  the  zeroes  of  C(z)  will  move  from 
t  the  origin,  towards  the  zeroes  of  A(z),  as  indicated  in  Figure  9.  The 

exact  trajectory  of  this  motion  can  be  plotted  using  classical  root  locus 
techniques  [7].  Note  that  the  zeroes  of  A(z)  are  shown  in  Figure  9  to 
be  on  or  very  close  to  the  unit  circle.  This  is  to  be  expected  for  narrow- 
band  line  spectra  and  for  pure  sine  waves. 

Several  conclusions  can  be  drawn  from  the  discussion  above:  (i)  For 
high  SNR,  no  pre-filtering  is  needed  since  C(z)  =  a  positive  constant; 

(ii)  For  very  low  SNR  C(z)  has  zeroes  near  the  unit  circle  which  means 
that  1/C(z)  will  most  likely  not  be  positive  real!  Thus,  pre-filtering 

A  A 

is  needed.  We  many  choose  either  D(z)  =  C(z)  or  D(z)  =  A(z) ,  since 

/v 

C(z)  *  A(z) .  The  choice  D(z)  =  A(z)  is  usually  preferred  since  the  esti¬ 
mates  of  the  AR  coefficients  {a. }  converge  much  faster  than M A  coefficients 

A  1 

(c±}. 

Preliminary  tests  of  the  algorithms  essentially  confirmed  these  con¬ 
clusions.  However,  serious  difficulties  were  experienced  in  the  case  of 
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Z-PLANE 


Z|,  =  zeroes  of  A(z) 


Figure  9:  Root  Locus  of  the  Zeroes  of  C(z) 
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narrowband  signals  in  which  case  A (2)  has  zeroes  very  close  to  the 
unit  circle.  Filtering  by  A(z)  led  the  algorithm  to  diverge  even  for 
reasonably  good  SNR's.  Further  investigation  indicated  at  least  two 
possible  causes  for  this  phenomenon: 

(i)  The  filter  A(z)  is  often  unstable,  i.e.,  A(z)  has  poles 
outside  the  unit  circle.  The  reason  is  that  since  A(z) 
has  poles  very  near  to  the  unit  circle ,  relatively  small 
estimation  errors  are  sufficient  to  make  A(z)  unstable, 
and  cause  the  algorithm  to  "blow  up”. 

(ii)  The  assumption  that  C(z)  =  A(z)  and  therefore  that 

A 

(A(z)/C(z)  -  1/2)  is  positive  real  is  only  true  for  very 
low  SNR's.  At  moderate  SNR's,  C(z)  may  be  quite  different 
from  A(z)  as  (31)  and  Figure  9  clearly  indicate.  Thus,  it 
would  be  preferable  to  find  a  filter  D(z)  that  is  closer 
to  C(z) . 

A  solution  which  addresses  both  of  these  issues  is  the  following: 

let 

/\ 

D(z)  =  A(kz)  (34) 

where  k  is  some  constant  smaller  than  one.  The  zeroes  of  A(kz)  are 
obtained  from  the  zeroes  of  A(z)  by  shifting  along  radial  lines,  as 
indicated  in  Figure  10.  The  new  filter  is  implemented  by  setting 

=  k1^  (35) 

since 

A  A  1  A  0  0  A  H  “H 

A(kz)  =  1  +  a,  kz  +  a.k  z  +  ...  +  a  k  az  a  (36) 

1  i  n 

a 

A  A 

The  modified  filter  A(kz)  is  more  stable  than  A(z),  since  its  roots 
are  further  away  from  the  unit  circle.  Furthermore,  by  a  proper  choice  of 
k,  these  roots  can  be  brought  closer  to  the  roots  of  C(z),  as  indicated 
by  a  comparison  of  Figures  9  and  10. 
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The  introduction  of  this  modified  pre-filter  greatly  improved 
the  convergence  properties  of  the  algorithms  for  moderate  and  low  SNR. 


Finally  we  should  note  that  typically  the  RML2  algorithm  is  used 

with  n,  ■  0  and  n  *  n  =  twice  the  number  of  sine  waves  expected, 
b  a  c 

Setting  n^  *  0  is  necessary,  since  the  inputs  u(t)  are  not  observable 
by  the  algorithm.  The  algorithm  is  also  used  in  another  mode  with 
na  =  n^  *  0  when  performing  TDOA  estimation  for  pure  sine  waves  in 
noise,  as  will  be  discussed  later. 


Figure  10.  The  Root  Locus  for  A(kz) 


I 


4. 


PERFORMANCE  EVALUATION 


The  MTS  algorithm  was  coded  and  tested  to  evaluate  its  performance 
for  different  types  of  signals  and  different  signal-to-noise  ratios.  The 
tests  so  far  have  been  restricted  to  a  single  fixed  target.  Two  aspects 
of  the  algorithm  were  studied  in  these  tests:  estimation  of  target  spec¬ 
trum  and  TDOA  estimation.  In  this  section,  we  present  some  preliminary 
results  which  indicate  the  type  of  performance  achievable  by  the  MTS 
algorithm.  It  should  be  emphasized,  however,  that  these  results  are  not 
conclusive;  more  testing  would  be  needed  to  establish  performance  bounds. 


4.1  Spectral  Estimation 


The  signals  used  in  our  spectral  estimation  experiments  were  sine 
waves  in  noise,  i.e.. 


m 

y(t)  =  "V  A  sin(2TTt/N . )  +  v(t) 
i=l  1 


where 

A.  =  amplitude 
N^  =  period 

v(t)  *  white  gaussian  noise 


(37) 


The  RML2  algorithm  was  used  to  identify  the  {  a^}  parameters  of  the 

received  signal  y(t).  The  final  estimates  a^  of  the  parameters  are 

then  used  to  generate  a  spectral  estimate.  In  our  simulation,  this  was 

done  by  generating  the  impulse  response  of  the  AR  model  1/A(z),  where 
~  ~  -1  ~  _na 

A(z)  =l+a^z  +  • * ’  +  an  z  ,  and  computing  its  power  spectrum. 

a 

Some  typical  results  for  two  test  cases  are  shown  in  Figures  11-20. 


Test  Case  #1 


The  signal  was  a  single  sine  wave  with  a  period  x  5.12. 

Assuming  that  the  MTS  algorithm  operates  on  a  1  Hz  bandwidth  (B  -  1  Hz) , 

fRECEQlfto  PASS 
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chis  corresponds  to  a  frequency  of  0.390  Hz.  It  should  be  remembered 

that  this  frequency  is  really  a  deviation  from  the  nominal  frequency 

used  in  the  bandwidth  selection  process  depicted  in  Figures  2  and  3. 

Thus,  we  actually  are  looking  at  an  expanded  picture  of  the  spectrum 

in  the  range  of  f  to  f  +B  Hz. 

oo 

The  algorithm  used  was  RML2  with  n  =  n  =*  2,  n  =*  0,  X(o)  =  .95, 

SC  D 

X^  =  .99.  In  all  cases  N=512  data  points  were  used.  This  corresponds 
to  an  integration  time  of  T  =  256  sec  =  4  minutes  (again  assuming  B=»l) . 

The  results  are  summarized  in  Table  1  and  Figures  11-15.  It  should 
be  pointed  out  that  for  the  low  SNR  cases  (-5  dB,  -lOdB)  the  algorithm 
really  requires  a  longer  integration  time.  However,  already  at  N=512 
points,  or  4  minutes  of  integration,  the  true  spectrum  starts  to  emerge. 

For  comparison  purposes,  we  have  included  in  the  figures  a  plot  of  a  con¬ 
ventional  FFT,  using  the  same  number  of  data  points.  Hanning  windowing 
was  used  where  indicated. 

Test  Case  #2 

The  signal  consisted  of  two  sine  waves  with  periods  N^=5.12,  N2=3.00, 
which  correspond  to  0.390  Hz  and  0.667  Hz.  The  same  algorithm  was  used  as 
in  Test  Case  ft  1.  The  results  are  summarized  in  Table  1  and  Figures  16-20. 
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True  (Windowed) 

-  Estimated  (Windowed) 

Figure  15A:  Estimated  and  True  Spectrum  --  SNR=-10dB 
P(o)=1,  KA=. 1 ,  IN IT=SPEC 


.  Estimated  (Windowed) 

-  True  (Windowed) 

Figure  19A:  Estimated  vs.  True  Spectrum  —  SNR=-5dB 
P(o)=.l,  KA=.5,  INIT=SPEC 


4.2  TDOA  Estimation 


Several  techniques  for  TDOA  estimation  based  on  the  MTS  algorithm 
were  implemented  and  tested.  The  first  and  most  straightforward  approach 
consists  of  estimating  the  coefficients  {b^}  of  an  MA  model  relating 
the  signals  y^(t),  y2(t)  at  the  output  of  two  receivers.  Let  us  assume 
that  the  signals  in  the  two  receivers  are  given  by 


y1(t)  =  xCt-D^  +  n^t) 


where 


y2(t)  =  x(t-D2)  +  n2(t) 


x(t)  ■  target  signal 


D^,D2  *  propagation  delays 

n^,n2  *  independent  measurement  noise  processes 
This  equation  can  be  rewritten  as 
y2(t)  -  yx(t-T)  +  n(t) 


where 


T  -  D2  -  Dx  -  TDOA 
n(t)  »  n2(t)  -  n1(t-t) 


Equation  (  39  )  represents  a  special  case  of  a  moving  average  model 


y2(t)  *  £  biy1(t-i)  +  n(t) 


with  b^  0  except  for  b^  •  1.  Thus,  estimating  the  model  parameters 
and  looking  for  the  largest  {b^}  will  indicate  the  value  of  the  TDOA. 
These  parameters  can  be  also  used  to  estimate  noninteger  values  of  the 
TDOA  by  a  proper  interpolation  technique,  as  discussed  in  Appendix  B. 
This  interpolation  technique  was  used  to  provide  estimates  in  two  test 


Case  //I 


A  second  order  AR  model  driven  by  white  noise,  with  a  spectrum 
given  by  Figure  21  .  The  algorithm  used  was  RLS  with  n  ■  n  *  0, 

3  C 

nfi  ■  7,  A(o)  «  .95,  \q  *  .99.  Some  typical  results  are  summarized  in 
Table  2.  The  true  value  of  the  TDOA  was  T  *  3.00. 


Case  #2 

A  fourth  order  AR  model  driven  by  white  noise  with  a  spectrum 
given  by  Figure  22.  The  same  algorithm  was  used  as  in  case  if  1.  The 
true  value  of  the  TDOA  was  3.00. 

TABLE  2 

TDOA  Estimation 
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A  somewhat  improved  version  of  this  approach  can  be  obtained  by 
using  the  MTS  algorithm  in  an  Adaptive  Line  Enhancer  (ALE)  mode  of  opera¬ 
tion.  The  algorithm  described  in  Section  3  provides  a  predicted  estimate 
of  the  input  signal  (see  Eq.  (14d)), 

y ( t+1)  -  «>(t+l)T  9 (t)  .  (41) 

The  estimated  signal  y(t)  provides  a  cleaner,  i.e.,  less  noisy  version 
of  the  received  signal  y(t).  This  is  illustrated  in  Figures  24  and  25 

A 

which  compare  the  power  spectra  of  y  and  y  for  two  test  cases.  Note 
the  significant  decrease  in  the  noise  levels  in  24B  and  25B. 

A 

Thus,  the  "enhanced"  signal  y(t)  can  be  used  as  an  input  to  the 
TDOA  estimation  algorithm,  as  depicted  in  Figure  23.  Initial  results 
have  indicated  some  improvement  when  this  method  was  used,  however,  more 
testing  is  necessary  before  final  conclusions  can  be  drawn 

A  second  approach  to  TDOA  estimation  is  based  on  "whitening"  the 
sensor  signals  using  the  MTS  algorithm  and  then  cross-correlating  to 
obtain  the  TDOA  estimate.  The  signal  whitening  is  achieved  by  using  the 
RML2  and  obtaining  the  residual  sequence  e  corresponding  to  the  input 
signal  yfc  (See  Figure  26. )  Correlating  the  residuals  gives  a  sharp  well- 
defined  peak  which  provides  a  better  indication  of  the  TDOA.  Some  typical 
examples  are  given  in  Figures  27A,  27B,  which  compare  the  correlation  func¬ 
tion  of  the  residuals  with  that  of  a  cleaner  version  of  the  data  obtained 

A  A 

by  using  the  predicted  signals  y^^*  This  approach  has  significant 
similarities  to  the  coherence  techniques  now  widely  employed  for  target 
detection  and  localization  [  9]. 
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Figure  25A:  Power  Spectrum  of  Two  Sine  Waves  in  Noise 
SNR=OdB 
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Figure  25B:  Power  Spectrum  of  Estimated  Signal 
RML2,  N=2048 
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— .  Residuals 

-  Predicted  Signal 

Figure  27A:  Correlation  of  Residuals  and  Predictions 
Signal:  a-j=-.606,  a2=.81,  SNR=10dB 

Algorithm:  RML2,  na=nc=2,  nB=0,  P(o)=lO,  KA=.8 

N=512,  \(o)= . 95 ,  A  =  99,  INIT-SPEC 


Residuals 


-  Predicted  Signal 

Figure  27B:  Correlation  of  Residuals  and  Predictions 
Signal:  a^.606,  a2=.81,  SNR=-dB 

Algorithm:  RML2,  na=nc=2,  Nb=0,  P(o)=l,  KA=.8, 

N=512,  A(o)=.95,  A  =.99,  INIT=SPEC 


A  third  approach  to  TDOA  estimation  is  based  on  the  idea  that  the 
residuals  e^(t)  provide  an  estimate  of  the  white  driving  process  u(t). 
Therefore,  it  is  possible  to  use  the  residual  £^(t)  computed  for 
sensor  It  1  together  with  the  received  data  (t)  in  sensor  It 2  as  the 
"known"  input  and  output  of  an  ARMA  model,  and  thus  apply  the  RLS  esti¬ 
mation  algorithm  to  find  its  parameters  (see  Figure  28). 


Figure  28.  TDOA  Estimation  by  Estimating  The  Input 
to  the  Spectral  Model 


The  residual  process  e^t)  is  in  fact  a  noisy  estimate  of  the 
input  process  u(t).  It  has  two  components:  one  due  to  the  measurement 
noise,  and  the  other  due  to  the  unpredictable  part  of  tne  signal.  In 
wideband  signals,  the  second  component  is  significant  and  we  may  expect 
e^Ct)  to  provide  a  reasonable  estimate  of  u(t).  However,  in  narrow- 
band  signals,  which  are  highly  predictable,  the  second  component  is 
small  and  e^(t)  is  a  very  noisy  estimate.  (In  fact,  for  pure  sine  waves 
the  second  component  vanishes ! ) 

These  statements  are  substantiated  both  by  theory  and  by  tests. 

We  found  that  for  pure  sine  waves,  the  residuals  eventually  converge  to 
the  measurement  noise,  and  no  longer  contain  information  about  the  signal. 
For  AR  processes  which  are  not  pure  sine  waves,  the  method  described 
above  worked  satisfactorily  in  sufficient  high  SNR.  The  more  narrowband 
the  signal,  the  worse  the  performance  obtained  for  a  given  SNR. 

The  last  approach  that  was  considered  for  TDOA  estimation  was  to  perform 
multichannel  (single  input-multiple  outputs)  parameter  estimation  using  an 
extension  of  the  RML2  algorithm.  One  form  of  the  multichannel  algorithm. 
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suitable  for  the  no  noise  case  (SNR*^) ,  was  implemented  and  tested  success¬ 
fully.  As  expected,  no  problems  occurred  in  the  no-noise  case.  The 
algorithm  requires  some  modifications  before  it  can  be  used  on  noisy  data 
(see  [ 2 1 ) . 

Computational  Requirements 

The  computational  requirements  of  the  MTS  algorithm,  as  any  other 
algorithm,  are  difficult  to  estimate  since  they  depend  strongly  on  a 
particular  implementation.  Furthermore,  a  major  part  of  the  computational 
load  is  due  to  data  handling,  I/O,  and  the  interactive  nature  of  our 
current  program.  However,  a  useful  indicator  of  the  amount  of  computa¬ 
tion  involved  is  given  by  counting  the  number  of  operations  (multiplies 

and  adds)  needed  to  compute  equations  (14)  and  (22) ,  which  constitute  the 

2 

basic  RML2  algorithm.  An  approximate  count  gives  ~(4n  +  5n  )  multiplies 
(where  n  =  the  number  of  estimated  parameters)  and  a  comparable  number  of 
adds,  per  single  update.  If  the  algorithm  operates  on  M  sensors  for  N 
data  points,  the  total  count  becomes 

No.  of  operations  ~(4n  +  5n2)MN  =  2(4n  +  5n2)MBT 

Assuming  a  typical  set  of  parameters: 

n  =  20,  M  ■  5,  B  =  10  Hz,  T  ■  1  sec. 

we  get  2x10^  operations  per  second.  It  should  be  emphasized  that  this 
figure  is  a  very  rough  estimate.  Alternative  forms  of  these  algorithms 
are  currently  available  which  are  more  efficient  (the  so-called  "fast" 
algorithms)*,  however,  they  were  not  implemented  at  this  stage  of  the 
development . 
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5. 


WORK  IN  PROGRESS 


As  mentioned  earlier,  the  results  presented  in  this  report  are  only 
preliminary.  We  are  continuing  our  investigation  in  two  principal 
directions : 

(i)  Algorithm  Development/Refinement 

The  experience  gained  in  testing  the  MTS  algorithm  leads 
us  to  believe  that  the  performance  achieved  so  far  can  be 
further  improved.  Some  of  the  specific  issues  which  are 
currently  addressed  include: 

•  improved  convergence  by  monitoring  the  stability  of 
the  filter  C(z),  and  adjusting  the  parameter  vector 
9  so  that  the  roots  of  C(z)  will  stay  inside  the 
unit  circle.  The  results  of  some  initial  tests  are 
depicted  in  Figures  29A-29D.  Note  the  very  substan¬ 
tial  improvement  that  was  obtained  compared  to  Figures 
15A,  19A,  20A,  and  the  fact  that  Figure  29D  corresponds 
to  SNR  =-15dB ! 

•  development  of  algorithms  that  incorporate  structural 
constraints  of  the  estimated  parameters  (e.g.,  the 
fact  that  the  {c^}  parameters  are  related  to  the  {a^} 
parameters  via  equation  (31)). 

(ii)  Algorithm  Testing  and  Performance  Evaluation 

After  developing  the  core  MTS  program,  we  are  now  in  a 
position  to  perform  a  more  comprehensive  set  of  tests 
to  study  the  performance  of  our  algorithms.  Specific 
issues  which  are  being  investigated  include: 

•  test  the  tracking  capability  of  the  MTS  algorithm 

on  synthetic  data  with  time  varying  target  parameters 
(TDOA  and  spectrum) . 
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Figure  29C:  Estimated  and  True  Spectrum  - 
N=1024,  Stability  Monitoring 
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Figure  29D:  Estimated  and  True  Spectrum  —  SNR=-15dB, 
N=2048,  Stability  Monitoring 


•  Test  algorithm  performance  under  a  variety  of 

conditions  including  multipath,  and  more  realistic 
(but  still  synthetic)  data. 

In  addition  to  this  work  which  is  part  of  Phase  I  of  the  project, 
we  are  also  studying  some  of  the  problems  to  be  addressed  in  Phase  II, 
i.e.,  the  extension  to  the  multiple  target  case.  This  extension  will 
involve  fitting  a  multi-input,  multi-output  (MIMO)  ARMA  model  to  the 
observed  data,  as  depicted  in  Figure  30. 

We  are  currently  studying  some  of  the  basic  problems  involved  in 
estimating  the  parameters  of  MIMO  systems  and  evaluating  the  modification 
required  to  adapt  our  current  MTS  algorithm  to  the  multitarget  case. 

Our  approach  to  the  multitarget  case  will  consist  of  two  steps,  as 
mentioned  in  the  introduction.  First,  we  plan  to  treat  the  no-noise  case. 
Some  of  the  fundamental  issues  that  need  to  be  addressed  are: 

•  Develop  an  algorithm  for  identifying  multi-input  multi¬ 
output  systems  with  unknown  inputs.  Current  techniques 
are  available  only  for  the  known  input  case.  Some  pre¬ 
liminary  work  was  already  performed  in  the  current  phase 
and  we  do  not  anticipate  any  major  difficulties. 

•  Investigate  the  special  structural  properties  of  the  MIMO 
case  (e.g.,  going  from  Left  Matrix  Fraction  Description  to 
Right  Matrix  Fraction  Description,  while  preserving  the 
structure -(see  Appendix  A,  Equation  (18)). 

•  Study  questions  of  identif iability  and  uniqueness  of  the 
MIMO  ARMA  model  and  their  relationships  to  achievable 
resolutions  (e.g.,  separation  of  closely  spaced  targets) 
and  to  the  discrimination  capability  of  the  MTS  algorithm. 
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SPECTRAL  MODEL 


PROPAGATION  MODEL 


X(z)  =  A-1 (z)  U { z ) 


A(z) 


A^Z)  0 
0  A2(z) 


X(z),  U(z)  are  Nxl  vectors, 


N  =  number  of  targets 


y(z)  =  B(z)  X ( z )  =  B(z)  A-1 (z)  U ( z ) 


H(z) 


y(z)  is  an  Mxl  vector,  M  =  number  of  sensors 


Figure  30:  Model  for  the  Multi  target  Data 


•  Implement  and  test  a  candidate  algorithm  with  emphasis  on  its 
tracking  ability.  The  objective  will  be  to  demonstrate  that 
after  track  initiation,  the  MTS  algorithm  can  provide  consistent 
tracks  of  several  targets. 

In  the  second  part  of  our  investigation,  we  will  extend  the  MTS 
algorithm  to  the  noisy  data  case.  Some  of  the  basic  issues  here  are: 

•  How  to  do  proper  prefiltering  for  the  MIMO  RML2  algorithm. 

•  Develop  the  positive  real  conditions  for  convergence  of 

the  MIMO  algorithm  and  find  a  way  of  improving  its  convergence 
(as  we  did  in  the  single  target  case). 

•  Implement  and  test  a  candidate  algorithm.  Run  a  variety  of 
test  cases  at  different  SNR’s  to  study  convergence  behavior. 

•  Use  the  experience  gained  to  develop  a  final  version  of  the 
MTS  algorithm  and  thoroughly  test  its  tracking  capability. 

This  second  step  will  probably  be  more  difficult  than  the  first,  and 
require  more  preparation  In  terms  of  developing  some  new  theoretical  results. 
However,  our  experience  from  the  first  phase  of  the  project  provided  us  with 
a  clear  understanding  of  the  difficulties  involved  and  we  feel  that  the 
goals  of  the  project  can  and  will  be  successfully  achieved. 

The  results  of  the  second  phase  of  the  MTS  project  will  provide  a 
significant  contribution  not  only  to  multitarget  tracking  but  also  to 
other  areas  of  interest  to  the  Navy  such  as:  adaptive  processing  of  multi¬ 
channel  signals  (noise  canceling,  adaptive  deconvolution,  adaptive  line 
enhancement,  etc.)  and  the  modeling  of  vector  time-series. 
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APPENDIX  A 


System  Identification  for  Multitarget  Tracking* 


*Published  in  the  Proceedings  of  the  13th  Asilomar  Conference  on  Circuits, 
Systems  and  Computers,  Pacific  Grove,  California,  November  1979,  and  was 
also  presented  to  the  National  Academy  of  Sciences  Panel  on  Applied 
Mathematics  Research  Alternatives  for  the  U.S.  Navy,  Washington,  D.C., 
November  2,  1979. 
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The  dae  saapLid  version  at  cr.ese  cutouts  will  ae  wtic- 


y.  ,'<)  *  :«(’<  -  3. .  -  c.  '<) 
/,[<)  «  oxC*  -  3.)  -  n,.<>  , 


:  -  sulc.  T.  *  Gjic,  T,  *  3._c  . 

Hoce  chac  Cha  delays  *, ,  T,  ara  issumea  co  be  inceger 

nuidples  of  cr.a  sampling  period,  ho  iiiiinulciss  arise 
wnan  c ^ ,  r,  are  non-vacagar  multiples  proviaen  c.nac 

cha  sampling  period  .c  is  proparly  cnosan.  Thus  point 
will  be  discussed  m  note  lecail  lacer. 

An  assumption  chac  is  ocean  naae  in  cna  coneaxt  of 
spectral  estimation  ;e.g..  cha  Maxi aua  Entropy  .lac.ccd 
[l-*i-[16i)  is  chat  cna  rictivac  signal  prccasa  xic> 
can  ba  raprasanced  as  an  aucoregrtssive  orocass  of  ociar 
i.  i.a. , 

n 

x(M  •  -£a.:u<c-i)  -  u(hj  .  ;2) 

i-L  * 

where  u(<0  is  a  whica  noise  driving  process.  This 
assumpeion  is  hoc  assanciai  co  our  approach  and  is 
incroducad  nere  for  siaplidey.  [Va  will  snow  lacer 
now  co  handle  sore  general  emlccers,  namely,  chose  with 
racionai  spectra.) 


£>—  y.Cc) 


y  c ) 


Figure  1.  Two  Seniors  ana  3ne  Target 
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Taxing  cm  a  :-tniii:jrj  of  Iq .  (2)  wa  gee 
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,d> 


•;<( :  I 

i  -M.  :  1  •  z*2'- 
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coefficient  of  o„,d) 

13 

;  2  )  :  *22 
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this  cass  is  1-1  •  3 . 
ocher  coefficients  art 
of  measurement  .moise. 
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Figure  2.  ::oce  that  :t>a  cirsc  significant  aon-raro 
coefficient  if  b.'c)  is  b.  .  ,  snc  che  fine  non-zero 

thus .  tne  TXa  in 
Li  non- zero  /aiues  if  tne 


vnere 
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Vritten  in  /actor  fora  the  tranafar  function  fraa 
che  driving  process  a  to  tna  sensor  outputs  y,  ,  y, 
is 


Td)  •  3(  z) X(z) 


-■nara 


3(2.' 


*(Z)  •  3(2) 


a(t) 


•J(r)  -X(z) 


v  4a) 


.*3’. 


1 


.  4b) 


;!ota  tnac  :.ne  cumarscor  of  cnas  tranafar  function  ton- 
tains  me  information  about  tr.a  cargac  location  i.a.  . 
tna  730A) ,  v<uie  tna  tanominacor  toot a ana  tna  ivnamics 
of  tba  tar gat  emitter  signal. 


:foce  tnac  it .  .  tan  be  rewritten  in  a  slightly 
tiff  tract  tors.  If  we  multiply  through  by  nt)  and 
cute  tha  inverse  r-ersnafora,  wa  gat 


v(k)  •  -  V  a,y('a-i)  * 
i-1  * 


3,u(k-l)  raCal  . 

(7a) 


iaaarica 

.1)  It  anould  oe  notch  that  tha  numerator  poly¬ 
nomial  found  by  such  nodal  fitting  oil!  not  oa  unique, 
tinea  vlchouc  naving  direct  measurements  of  :c;  <c)  tha 
absolute  dalays  .i.a.,  tha  dagraas  2.,  2,  of  d.na 

polynomials  b.  ,  o,,'  cannot  oa  datarmlnad.  However, 
tha  dlffarenca  in  the  dagraas  of  tha  polynomials  o., 
b,  will  be  unique  ana  equal  to  tha  730A,  1.  ,*3.-3, , 

wnlcn  proviqas  tna  t as it a a  information  aoout  the  aourqa 
bearing. 

(2)  Figure  2  also  proviqas  an  indication  of  vnac 
will  happen  if  tha  dalays  art  non-integer  nuitiplas  of 
tha  sampling  period:  instead  of  naving  a  single  non- 
taro  coefficient  associated  with  a  given  delay,  wa 
will  have  two  large  coefficients  wnose  relative  nagni- 
tuaes  reflect  tow  close  tha  real  celay  is  to  the  delay 
represented  by  thac  coefficient .  'or  example,  if 

act  <  -,  <  (kvl)de,  we  nav  axpect  both  b.  .  and 
-  *  -  -  *  a 

b,  ....  to  ba  non-tars.)  evidently .  as  long  as  the 
I ,  i,  a-el) 

sampling  races  are  sufficiently  nign  compared  :s  che 
oandwidth  of  the  underlying  orocess,  the  same  aoproacn 
will  wore  for  noa-inceger  tiaa-deiays. 

(3)  tha  approach  described  hare  can  also  .nandie 
pulcipach.  In  tha  case  of  multipath  b. (z)  and 


wnara  ad)  is  a  correlated  noise  process  given  by 

a(’d)  •  a  (’a)  *  £  a.n(k-i)  .  (7b) 

1-1  * 

la  oc.ner  words,  the  ouepue  vector  y(k)  can  be  repre¬ 
sented  as  an  auto— regressive  moving- average  (ASHA) 
process  of  order  n. 


b,d))  will  have  several  large  toafficlancs  torras- 

ponding  to  the  direct  path  deiay  and  tha  multipath 
delays.  Since  the  direct  path  nas  tha  shortest  delay, 
tha  firsc  large  coefficient  of  b,(s)  will  correspond 

to  the  direct-path  delay.  Thus,  the  730A  can  be  easily 
evaluated  from  the  (3.(-  coefficients  even  in  the 

i 

presence  of  multipath! 


equation  (7)  immediately  suggests  a  method  for 
estimating  730 A:  find  a  sac  of  coefficients  ‘a, ,3, ; 

thac  best  fits  the  lace  yd)  in  the  seen  square 
error  sense.  7his  can  be  done  by  varloue  parameter 
estimation  techniques  to  be  ilscusaed  in  a  lactr  sac- 


tion.  Once  toe  estimates  (a^,S^> 

art  found, 

:h« 

TOOA  can  be  evaluated  by  looking  at 

3, . 

To  see  this  hore  clearly,  consider  the  axampls 

of  Figure  1,  with  3.  *1,  3,  *3,  !. 

(k)  •  x(k-i)  *  a.  (k)  , 

(3a) 

y,d)  ■  x(*-*)  *  a,(k)  , 

(3b) 

which  means  720 A  •  3 .  In  this  casa , 
Ic.  (4), 

according 

to 

b. (2)  ■  2**  , 

(9a) 

.  •<* 

3,x2)  •  C2 

(9b) 

7hus,  if  wa  estimate  the  coefficients  of  b,(z),  b,(r) 

by  perfornihg  sn  aSHA  fitting  on  the  data  sac  (yd);, 
wa  will  expect  to  see  the  situation  depicted  in 


THE  XDiri7.UtCE7  CASE 

The  ASHA  modeling  approach  can  ba  eaeily  extended 
to  the  multitarget  case,  Here  che  system  consisting 
of  eargses  and  sensors  will  be  represented  by  a  multi- 
input,  multi-output  transfer  function  (i.a.,  ASHA 
model).  A  simple  example  is  depicted  in  Figure  3. 

Tha  equation  describing  chs  vector  of  measured 
deta  yd)  is  given  by 


1 


Figure  2.  Numerator  Coafflclants 
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one  re 


•  3d  )\izi 


:.o*) 


3c!  •  .3,  z)  ,3„,c)  i  ■  (10b) 


If.  u  feiore.  ve  assume  mac  the  signals  emitted  by 
:ic;icj  zi  ir*  sucocegrtssive  processes,  then 


fa.  ’z) 

, 

’  i.  (2)  ' 

A( Zi  *  1 

,  w( 2)  - 

- 

l  3 

^(2) 

thus,  me  output  oi  ".Si  sensors  in  the  nuicicargec 
rasa  ran  oe  writ c  an  n 

•'(J)  •  3(1!  A*1-!!!  3(1)  -  M< z)  .  (12) 


the  ruascion  oi  finding  the  target  “OA’s  '.1  equivalent 
: o  oroolea  si  escmaciag  me  roeffi  atncj  .3..  ,A,  } 

of  one  transfer  function  3(:)A  '  ’z)  ■  2n  ocher  words, 

rind  i  aodii  3  ( c )  A  4(2)  that  sill  rase  fit  the  avail- 
lol*  iaca  tnce  the  nodel  has  scan  round,  che 

location  o i  me  target  sill  p a  round  Oy  examining  me 
roaiflcienca  in  the  appropriate  column  oi  3(:)  .  and 
me  cargec  ipeccra  ran  Pa  round  from  the  ippropriaca 
alaaanea  or  A(zj. 


Actually,  ip .  (12)  la  not  quite  me  fora  sa  gat 
scan  Aitf!A  aooeliag  la  performed.  since  me  auiti-inpuc. 
auiti-ouepue  ARMA  nodal  Pas  me  fora 


y(’*;  *  - ^ a, y(k-t;  -  t3,u<k-i)  . 

i-l  *  i-l  • 

(13) 

ig  2-transtoraa  will  giva 

A(Z)  ?(S)  •  3(2)  3(2)  . 

(14) 

?(2)  •  A‘l(2)  3(2)  2(2)  , 

C-5) 

~(Z)  •  I  “  T  A, 2**  , 

i-l  ' 

'  15) 

3(2)  *  Tl,  c"1  . 

(It) 

1-1 


3y  comparing  .  11.  ana  10/  sa  rocs  mac 
rinanca  A,  .3,  are  roc  praciaai v  rhcaa  oi 


rsar- 

trans- 


:'t.-  function  3C;  a”*C*j  .  in  fact,  may  ire  relacaa  cy¬ 
me  apuacion 

A  3(r;  •  3(2)  a  ~\Zi  13/ 


Thus ,  If  sa  use  some  parameter  estimation  rscnnipua  rs 
fit  in  AfUiA  nodal  A,  ,3,  to  me  rata  yik)  ,  sa  sill 

-tave  to  perform  aitarsardi  me  acap  of  computing  A,, 3,. 
froa  !a,.3,;,  and  than  evaluate  che  ttOA’s.  Various 

esenniques  are  available  for  performing  mesa  rompuca- 
nona  '  11 1  -  ( 12  ] . 

!foce  mac  m  tpa  proposed  approach,  mera  la  no 
read  to  "laoal"  the  targets  at  each  imp.  It  is  neces¬ 
sary,  oi  toursa,_:o  escaoiisn  initially  me  laoeiing  of 
me  loloans  oi  3(zi  or  A(ti),  i.e.,  ram  mine 
smen  column  mall  rarer  to  snitn  target.  Attarsart. 
tPa  esc  mac  as  oi  the  niiMA  coefficients  aca  rpcatan  at 
aacn  r me  step  a  according  to  me  new  lata  sector 
yU)  • 


The  proposed  aporoacn  periocas  i  recursive  giooai 
estimation  process,  i.a.  .  for  ail  targets  over  me 
tPosen  tiaa  interval.  The  algorithm  aucsaacicaily 
tries  to  fit  a  sec  oi  coefficients  snitn  'explain''  in 
the  Peat  say  all  me  available  iaca.  3inca  mate 
coefficients  toncain  the  TtOA  mfomation  for  ail  me 
targata  /as  sail  as  their  spectra)  ,  se  get  a  sec  oi 
r  oasis  cant  esemaces  oi  all  targtc  locations,  Pa  sac  on 
all  the  availaole  iaca.  if  an  retinal  parameter  esci- 
saeor  13  used,  me  resulting  TtOA  estinatas  are  truly 
opernai. .  In  the  standard  tracking  approacnes  iiseussan 
in  rha  Introduction,  only  suDoociaal  parameter  esci- 
nates  can  oe  oocainec,  since  each  target  location  Is 
first  esc. maced  ladivicually  and  me  cracpriag  aigoriina 
then  actampcs  an  appropriate  labeling  oi  these  esti- 
nacas. 

finally,  it  snould  Pa  noted  mac  finding  che  esci- 
natas  oi  the  \  A^;  coefficients  is  equivalent  to  per¬ 
forming  sinuicanaous  soactral  estimation  oi  ail  me 
targets,  with  automatic  line  association,  the  iizsc 
is  true,  since  me  spaccrum  can  p«  computed  directly 
from  the  autoregressive  coaiticiancs  fuse  as  for  chs 
Maximum  Incropy  dec  hod  [  13  j  —  C  lj  i  -  the  latter  is  true 
for  the  same  reason  that  ao  relabeling  oi  the  70CA ' s 
is  requirsd. 

SSCOTSIVE  PARAdETZR  iSTIMAttON 

The  approach  outlined  in  the  previous  sections 
depends  on  our  ability  to  tompuce  rha  aRMA  coefficients 
given  a  sec  oi  neasureaents .  this  type  oi  problem  has 
Paen  widely  studied  in  the  general  loncaxc  oi  parameter 
estimation  and  in  tht  note  specific  context  oi  identi¬ 
fying  svseen  aodels  from  iaput/oucpuc  naasuranents  ; 15 i — 
[241. 


the  least  squares  paraaactr  astiaaclan  problem  is 
usually  formulacad  as  follows:  given  a  sec  oi  measure- 
nenta  {y(k)  ,u(«c)  ).em^  find  the  coefficients  {a,  ,3,; 

that  will  minimize  tha  naan  square  error  £,iy(iO-y(ic)  !  * 
wnara  k 

?0t>  *  *  ^  A,7(k-l)  *  ^3,u(i-i)  .  ( 19) 

i  *  i  * 

tha  vector  y(k)  is  tns  value  predicted  Pv  tha  ARPIA 
nodal,  for  tha  naasuraaant  at  time  s. 
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•  -  v  4 r £  ,  1..9  >w*uC*jn  * 

::.^3  pcoiLdra  .j  fairly  icrai^nciacvara  iaa  typira-  il ,. - 
rithms  -»n  :*  four-a  in  [lo;,  lir),  111),  !-I|.  The 
situation  Li  jomewnac  r.ort  tom? Licacas  men  . :  Li  no: 
possible  ::  measure  1;  i,  .  Sowever.  ay  assuming  :;,ic 
:.k,  Li  i  iacuencs  of  mcepencenc  '"units"'  random 
rariaoles  visa  l3.:  /anur.ca .  it  lj  jell*.  possible  to 
asciaace  tne  noaei  parameters.  This  case  ls  usuailv 
:i:ac;ta  to  is  "tne  :u<  of  torralacac  tesicuais,  "  anc 
several  Lscnniguas  *sv»  oeen  suggescea  for  its  aolu- 
;ioa.  -or  itciiii.  see  the  iu rvev  by  Ascrom  ' L5 i . 


Mors  rscsnciy  1  new  10 prose a  ass  ossa  aeveiopea 
oy  Mori  ,'13)-;20!  mica  provides  efficient  forms  or  the 
jo-called  exact  recursive  lease  jausrts  algoricnm*. 
Thess  isu  forms  asvs  tas  idcsa  aavsneages  or  computa¬ 
tional  efficiency  laa  fisc  peramecer  tracking  -.apaoii- 
iS7  .'231,  ^24,' .  Ths  last  proparry  is  laporrsac  for 
cracking  moving  sines  Ehs  ARMA  aodsi  iorres- 

ronding  :o  such  targets  has  Lias  v arying  parameters. 
Thtss  a-goricnas  srs  also  capaol*  of  nandling  r.ansca- 
tionar'  source  sna  toiss  processes. 
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APPENDIX  B 


TDOA  Estimation 
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Let 


y^t)  =  x(t)  +  n^(t) 
y2(t)  =  x(t+T)  +  n2(t) 


(Bl) 


-epresent  signals  received  by  two  different  sensors.  The  noise  processes 
n^(t),  n2(t)  are  assumed  to  be  white,  and  independent.  The  two  signals 
are  related  by 


y2(t)  =  yx(t+T)  +  n(t) 


(B2) 


where 


n(t)  =  n2(t)  -  n1(t+T). 

The  noise  process  n(t)  has  a  variance  equal  to  the  sum  of  the  variances 
of  n^  and  n2.  The  sampled  values  of  y^,y2>n  will  be  denoted  by 
y^kAT),  y2(kAT),  n(kAT).  Assuming  that  the  sampling  interval  AT  is 
adequately  small  for  x(t),  we  have 

+00 

y,  (t)  =  £  sinc(t-kAT)  (B3) 

k=-°° 


where 


sinc(t)  *  sin(Trt/AT)/(Trt/AT) . 


(B4) 


Let 


x  =  IAt  +  At  ,  o  <_At  <  At. 

+00 

y2(iAT)  =  Yj  yxCkAT)  sinc[(i+«.-k)AT  +  Ax]  +  n(iAT)  (B5) 

k“-<» 
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Without  loss  of  generality,  we  can  set  AT  »  1,  and  make  a  change 
of  variables  i-k  =  n,  which  will  give 

+00 

y2(i)  =  b^Ci-n)  +  n(i)  (B6) 

n=-°° 


where 


b  =  sinc(n  +  l  +  Ax)  (B7) 

n 

Thus,  the  time  series  y2(i)  is  related  to  y^(i)  by  a  moving  average 

(MA)  filter  with  coefficients  as  given  by  Equation  (B7) .  In  practice, 

we  will  consider  only  a  finite  number  (n  )  of  terms  in  the  sum  (B6) . 

b 

The  coefficients  bn  can  be  considered  as  the  samples  of  a  function 
sinc(n  +  l  +  Ax)  which  achieves  a  maximum  at  n  +  1  +  Ax  =  0.  Hence, 

A 

given  the  coefficients  b^,  the  delay  X  is  the  value  which  maximizes 
the  function- 

% 

b(x)  =  T]  b  sinc(x-n)  .  (B8) 

n-1  n 

In  our  experiments,  we  used  a  search  algorithm  to  find  the  value  x 
which  maximizes  b(x).  Some  typical  results  are  summarized  in  Table  2, 
Section  4.  A  similar  approach,  which  uses  a  different  type  of  estimation 
algorithm,  can  be  found  in  [8]. 
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The  MTS  algorithms  were  implemented  on  SCI's  VAX.  The  programs  are 
written  in  FORTRAN  and  are  fully  interactive.  Plotting  capabilities 
include  a  Tektronix  display  and  character  displays.  The  interactive 
program  allows  easy  changes  of  test  cases  (target  spectra,  signal-to- 
noise  ratio)  and  algorithm  parameters  (type  of  algorithm,  model  order), 
as  well  as  convenient  program  modification. 

The  following  pages  present  an  example  of  the  program  parameters 
under  our  control  and  some  typical  plots  obtained  for  a  sample  test  case. 
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